program mkhisq ! -- makes histograms from list-mode data file.
  !
  ! DESCRIPTION:
  !
  ! This is a speed-enhanced version of mkhis. Cumbersom, painful
  ! checks about validity of the list-mode data is reduced as much as
  ! possible. It only checks the pair-flags and rejects the event on
  ! the same plane.
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use kuty
  use suty
  implicit none
  interface
     subroutine decvq(pe,ne,pb,nb,pneo,pnbo)
       integer*8,intent(in),value::pe,ne,pb,nb,pneo,pnbo
     end subroutine decvq
  end interface
  integer,parameter::MAXPATHLEN=1024
  integer*8,parameter::fbs=5898240
  integer*8,parameter::nemax=fbs/4
  integer*8,parameter::nlmdbufmax=737280
  integer*4,allocatable,target::lmdbuf(:,:)
  integer*8 nlmdbuf_wrote, nbuf_read
  !
  integer,parameter::nd=8         ! total number of blocks
  integer,parameter::ndh=nd/2     ! number of blocks in a plane
  integer,parameter::mgap=4       ! number of crystals in a gap
  integer,parameter::nt=64        ! number of crystals in transaxial
  integer,parameter::na=20        ! number of crystals in axial
  integer,parameter::nx=na*ndh+mgap*(ndh-1) ! number of crystals in x
  integer,parameter::ny=nt        ! number of crystals in y
  integer,parameter::C2=ndh*(na+mgap)
  integer,parameter::nxy=nx*ny
  !
  integer*8 pbuf
  integer*4,allocatable::buf(:)
  integer*2,allocatable,target::hists(:,:,:)
  integer*2,pointer::hist(:,:),histp(:,:),histd(:,:)
  character*(MAXPATHLEN) fin,str,flist
  character*(MAXPATHLEN),target::fout,foutp,foutd
  character*(MAXPATHLEN),pointer::s
  integer istat,len
  integer*8 n,ie,fs,ne,nbuf,nb,nread
  integer*4,pointer::at1,at2,pd
  integer if,nf,ia,narg  
  integer ndelay,npair
  logical useflist

  narg=command_argument_count()
  if(narg<3) then
     call usage
     stop "*** too few arguments."
  end if
  
  fout=""
  foutp=""
  foutd=""
  flist=""
  nf=0
  nullify(s)

  do ia=narg-1,2,-2
     call get_command_argument(ia,str,len)
     if(len==2) then
        select case(str(2:2))
        case("o")
           s=>fout
        case("p")
           s=>foutp
        case("d")
           s=>foutd
        case default
           call usage
           stop "*** invalid option."
        end select
        if(associated(s)) then
           call get_command_argument(ia+1,s)
           nullify(s)
        end if
        nf=ia-1
     else
        if(fout/="".or.foutp/="".or.foutd/="") exit
        call usage
        stop "*** please specify outfile correctly."           
     end if
  end do

  call get_command_argument(1,str)
  if(str=="-f") then
     call get_command_argument(2,flist)
     nf=huge(1) ! <
     open(unit=3,file=flist,iostat=istat)
     if(istat/=0) stop "*** cannot open flist"
     useflist=.true.
  else
     useflist=.false.
  end if

  allocate(hists(0:nxy-1,0:nxy-1,0:1),stat=istat)
  if(istat/=0) stop "*** cannot allocate memory"
  allocate(lmdbuf(4,nlmdbufmax),stat=istat)
  if(istat/=0) stop "*** cannot allocate memory"
  histd=>hists(:,:,0)
  histp=>hists(:,:,1)
  hists=0
  ndelay=0
  npair=0
  nb=0
  nbuf=0
  pbuf=0

  if=0
  file: do
     if(useflist) then
        read(3,"(a)",iostat=istat) fin
        if(istat/=0) exit
        if(len_trim(fin)==0) cycle
     else
        if=if+1
        if(if>nf) exit
        call get_command_argument(if,fin,len)
     end if
     fin=trim(adjustl(fin))     
     write(*,"(a,x,$)") "Opening file: "//trim(fin)
     open(unit=1,file=fin,status="old",access="stream",iostat=istat)
     if(istat/=0) stop "*** cannot open file."
     fs=fsize(fin)
     write(*,*) trim(i2a(fs))//" byte"
     if(fs==0) cycle

     ne=fs/4
     n=min(ne,nemax)     
     if(nbuf<n) then
        if(nbuf>0) deallocate(buf)
        allocate(buf(n),stat=istat)
        if(istat/=0) stop "*** cannot allocate memory"
        nbuf=n
     end if

     nread=0
     read: do        
        read(1,iostat=istat) buf(1:n)
        if(istat/=0) stop "*** error during read."
        nread=nread+n
        nb=nb+n
        pbuf=loc(buf)

        dec: do
           call decvq(pbuf,n,loc(lmdbuf),nlmdbufmax,&
                loc(nbuf_read),loc(nlmdbuf_wrote))
           npair=npair+nlmdbuf_wrote
           sort: do ie=1,nlmdbuf_wrote
              at1=>lmdbuf(1,ie); at2=>lmdbuf(2,ie)
              pd=>lmdbuf(3,ie)
              if(pd==0) ndelay=ndelay+1
              hists(at1,at2,pd)=hists(at1,at2,pd)+1
           end do sort

           if(nbuf_read==n) exit
           pbuf=pbuf+4*nbuf_read
           n=n-nbuf_read

        end do dec

        if(nread==ne) exit          
        n=min(ne-nread,nemax)

     end do read

     close(1)

  end do file

  write(*,*) "Number of files read:    ",nf
  write(*,*) "Number of bytes read:    ",nb*4
  write(*,*) "Number of accept events: ",npair
  write(*,*) "Number of random events: ",ndelay
  
  if(foutp/="") call write_hist(histp,foutp)
  if(foutd/="") call write_hist(histd,foutd)

  hist=>histp
  write(*,*) "Max Prompt: ",maxval(histp)
  write(*,*) "Max Delayed:",maxval(histd)
  hist=histp-histd
  where(histp<0) histp=0
  write(*,*) "Max True:   ",maxval(hist)

  if(fout/="") call write_hist(hist,fout)
 
contains
  
  subroutine usage
    write(*,100)
100 format(x, "Usage: mkhis {infiles | -f flist}",/,&
         8x,"{-o outfile_true | -p outfile_prompt | -d outfile_delayed}")
  end subroutine usage
  
  subroutine write_hist(h,outfile)
    integer*2,intent(in)::h(:,:)
    character*(*),intent(in)::outfile
    integer istat
    character*128 s
    write(*,"(a,$)") "Writing "//trim(outfile)
    open(unit=1,file=outfile,access="stream",status="replace",iostat=istat)
    if(istat/=0) then
       write(*,*) "*** error during open"
    else
       write(1,iostat=istat) h
       close(1)
       if(istat/=0) then
          write(*,*) "*** error during write"
       else
          s=trim(i2a(nx))//"x"//trim(i2a(ny))
          write(*,*) "size: "//trim(s)//"x"//trim(s)//"x"//trim(i2a(sizeof(h(1,1))))//" = "&
               //trim(i2a(sizeof(h)))//" bytes"
       end if
    end if
  end subroutine write_hist

end program mkhisq
